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Abstract 

The accuracy of Green Function Monte Carlo (GFMC) simula- 
tions can be greatly improved by a clever choice of the approximate 
ground state wave function that controls configuration sampling. 
This trial wave function typically depends on many free parame- 
ters whose fixing is a non trivial task. Here, we discuss a general 
purpose adaptive algorithm for their non- linear optimization. 

As a non trivial application we test the method on the two 
dimensional Wess-Zumino model, a relativistically invariant super- 
symmetric field theory with interacting bosonic and fermionic de- 
grees of freedom. 



1 Introduction 



The traditional algorithms for numerical simulations of Lattice Field The- 
ories are based on the Lagrangian formulation [Q] that follows Feynman's 
idea of representing quantum amplitudes in terms of sums over classical 
paths. An interesting alternative is the Hamiltonian framework of Kogut- 
Susskind J2J focusing on the Hamiltonian as the generator of the temporal 
evolution. To improve numerical efficiency in the evaluation of physical 
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observables, analytical approximations of the ground state wave function 
can be successfully exploited. This well known procedure is called Impor- 
tance Sampling || and its key problem is to determine the optimal trial 
wave function (TWF) within a given class depending parametrically on a 
set of free parameters. 

In this contribution, we review a recently proposed algorithm that 
solves this problem. It is built upon a standard Green Function Monte 
Carlo (GFMC) simulation algorithm and is based on a non trivial feed- 
back between Monte Carlo evolution and TWF modifications. 

To discuss the method, we present some novel results concerning the 
calculation of the ground state energy of the supersymmetric N — 1 Wess- 
Zumino model in 1 + 1 dimensions. 

2 Review of Green Function Monte Carlo 

GFMC algorithms can be regarded as numerical implementations of the 
Feynman-Kac formula |4j] . In the simple case of Quantum Mechanics, this 
formula claims that, given the potential V(q) and an initial wave function 
ipo(q), the function 

^(q, t) = E L Q (q + W t ) exp - J* V(q + W s )ds^ , 

(where W t is the Wiener process) solves the Euclidean Schrodinger equa- 
tion 

d t ij(q,t) = ^Aij(q,t)-V(q)?P(q,t) 

with initial data ip(q, 0) = ipo(q). 

A translation of this formula into a numerical algorithm is straight- 
forward and quantum expectation values over the ground state can be 
expressed as statistical averages over an ensemble of weighted walkers 0. 
As is well known, the weight variance of the ensemble members must be 
controlled in some way. Here we adopt the Stochastic Reconfiguration 
algorithm with fixed population size 0. The finite population size bias 
as well as walker correlations vanish with increasing number of walkers. 

Weight fluctuations are closely related to the noise of measurements 
and are due to the fact that V is not constant along walker paths. To 
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significantly improve accuracy one needs to reduce these fluctuations. Im- 
portance Sampling is a common strategy to achieve such a noise reduction: 
the original Hamiltonian H = \p 2 + V(q) is transformed into 

H = e F He~ F = ^p 2 +ip-VF + V, 

V = V - \&F- ^(VF) 2 . 

H is not canonical, but still allows simulations with minor modifications 
with respect to the F = case [EL The diffusion of the walkers is driven 
by a drift term proportional to V-F while the potential V(q) is replaced 
by a new potential V that depends on the trial wave function through F. 
Roughly speaking, a true improvement is reached when V is more constant 
than V along random paths in some way we are going to discuss later on. 



3 Adaptive Optimization of the Wave Function 

In this Section, we show how the trial wave function F can be optimized 
dynamically within Monte Carlo evolution. To this aim, we consider a 
parametric TWF exp -F(q, a) depending on the state degrees of freedom q 
and on a set of free parameters a. After iV Monte Carlo steps, a simulation 
with a population of K walkers furnishes a biased estimator Eq(N, K, a) 
of the ground state energy E that we shall take as a representative ob- 
servable. Eq is thus a random variable such that 



(E (N, K, a )) = E + ^ + (K- a ), a > 



and 

Var E (N, K, a) 



c 2 (K, a) 



N 



where (•) is the average with respect to Monte Carlo realizations. 

In the K — > oo limit, (Eq) is therefore exact and independent on a. The 
constant C2(K, a) is related to the fluctuations of the effective potential 
V and is strongly dependent on a. The problem of finding the optimal F 
can be translated in the minimization of Oz(K, a) with respect to a. 
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The algorithm we propose || performs this task by interlacing a Stochas- 
tic Gradient steepest descent with the Monte Carlo evolution of the walk- 
ers ensemble. At each Monte Carlo step, we update a n — > a n+ i according 
to the simple law 

a n+ i = a n - r] n V a V&Y £n V 

where £ n is the ensemble at step n and {r/} is a suitable sequence, asymp- 
totically vanishing. 

A non-linear feedback is thus established between the TWF and the 
evolution of the walkers. The convergence of the method can not be easily 
investigated by analytical methods and explicit numerical simulations are 
required to understand the robustness of the algorithm. In reference f|, 
examples of applications of the method with purely bosonic or fermionic 
degrees of freedom can be found. In the next Section, we shall address a 
model with both kinds of fields linked together by an exact lattice super- 
symmetry. 



4 The N = 1 Wess-Zumino Model in 1 + 1 Dimensions 
and Supersymmetry Breaking 

The model we consider is a lattice version of the N = 1 Wess-Zumino 
model previously studied in 0. On each site of a spatial lattice with 
L sites, we define a real scalar field {<p n } together with its conjugate 
momentum {p n } such that [p n , <p m ] = —i5 n)m . The associated fermion 
is a Majorana fermion {ip a ,n} with a = 1,2 and {ip a ,m V^m} = 8 a ,b8n,m , 
4>l n = ip a ,n- The fermionic charge 



Q =Y1 Pn01,n - ( V? " +1 „ — 1 + V((p n )j 1p2,n 
n=l L V ^ / 



with arbitrary real prepotential V((p), can be used to define a semi-positive 
definite Hamiltonian H = Q 2 . Its positive eigenstates are automatically 
paired in doublets connected by Q; only a Q-symmetric ground state with 
zero energy can be unpaired. H describes an interacting model, symmetric 
with respect to Q; its ground state breaks the symmetry if and only if its 
energy is positive. We stress that spontaneous supersymmetry breaking 
can occur even for finite L, because tunneling among degenerate vacua 
connected by Q is forbidden by fermion number conservation. 
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We can replace the two Majorana fermion operators with a single Dirac 
operator satisfying {Xn,Xm} = 0, {Xn,xln} = 8 n,m 0- The Hamiltonian 
takes then the form 



and describes canonical bosonic and fermionic fields with standard kinetic 
energies and a Yukawa coupling. 

The theoretical problem that we want to address is dynamical breaking 
of supersymmetry in the above kind of models. At tree level, supersymme- 
try breaking occurs if and only if V(ip) has no zeroes. In 1 + 1 dimensions, 
this picture is known to be unstable with respect to radiative corrections 
as shown by the analysis of the one-loop effective potential of the scalar 
field (p. From the non-perturbative point of view, the ground state energy 
E is the crucial quantity because it vanishes if and only if we are in the 
symmetric phase. The choice of a Hamiltonian Monte Carlo algorithm 
is then natural since E is the simplest observable that can be measured 
with such a technique. 

We study two specific examples: the cubic prepotential V((p) = (p 3 , 
where unbroken supersymmetry is expected, and the family of quadratic 
prepotentials V((p) = X + ip 2 . In the latter case, at tree level supersymme- 
try is broken for Ao > 0. Perturbative calculations and previous numerical 
results [0] suggest that this conclusion does not hold at the quantum level 
and provide evidence that supersymmetry is broken in the lattice model 
for Ao greater than some negative Aq < and is recovered for Ao < Aq. 
On the other hand, the strong-coupling limit suggests E > for all Ao, 
decreasing exponentially for Ao — > — oo ||; moreover, one can show that, 
at fixed finite L, E is an analytical function of Ao; this would imply that 
a symmetric phase is only possible in the L —>■ oo limit. 

4.1 Simulation Algorithm: Brief Description 

The Monte Carlo evolution in GFMC algorithms approximates the Hamil- 
tonian evolution semigroup {e~ tH } t >o- The bosonic and fermionic degrees 




Hb(p,<p) + h f (x,x,v) 
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of freedom can be split by writing 

e -W = e -l/2 (3H Be -(3H F e -l/2 /3H B + 0(f3 3 ), 

and the /3 — ^ extrapolation is numerically computed at the end. We 
implement the evolution associated to Hb to second order in (3 @ and the 
evolution associated to H F exactly [|1(| . 



About the choice of the trial wave function, we propose the factorized 
form 

where \^o)b ® |^o)f is the ground state for the free model and 

Sb = EE* ^ = EC- 1 )" (xlx» - 5) f>K- 

n fc=l n v z/ fc=i 

The degrees <ie and e2p must be chosen carefully to achieve convergence 
of the adaptive determination of {a B ,a F }. 

In the regime of weak coupling perturbation theory, the choice of peri- 
odic boundary conditions does not break supersymmetry when Lmod4 = 
0. Under this condition, there is an even number of fermions, L/2, in the 
ground state. A sign-problem then arises due to boundary crossings since 
they involve an odd number of fermion exchanges. To avoid such a diffi- 
culty, we shall adopt open boundary conditions. With this choice, L needs 
just to be even to assure a supersymmetric weak coupling ground state. 
In the following, we shall restrict ourselves to the case Lmod4 = 2 and 
to the sector with L/2 fermions that contains a non-degenerate ground 
state, with zero energy at all orders in a weak coupling expansion. 



4.2 Numerical Results 

We choose ds — dp — 4; due to the symmetries of the model with 
quadratic V((p) = Ao + tp 2 , we impose parity constraints on the trial 
wave function, setting af , a^, af, af to zero (af , , af , af in the 
cubic case), and and determine by the adaptive algorithm the remaining 
four free parameters. We work on a lattice with 10 spatial sites and vary 
both K and A . In principle, we also extrapolate any result to the (3 — > 
limit. However, the results we present are all obtained with a fixed value 
(3 = 0.01 since we checked that a reduction of this parameter by a factor 
2 does not change appreciably the results. 
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Figure 1: Adaptive Optimization of the four coefficients appearing in the trial 
wave function for the quadratic model. The constant Ao is 0.0 and the ensemble 
size is K = 100. 



Beginning with a small ensemble of K — 100 walkers, we determine 
adaptively the four coefficients {af , af , af , af} starting from zero values. 
A typical run is shown in Fig. ([]]). For larger K, we start from the values 
of a obtained in a run at the nearest smaller K. A summary plot of 
{af , af , af , } as functions of K and Ao is shown in Fig. (|2|). 

With the optimized TWF parameters, we determine the energy con- 
tributions coming from the bosonic and fermionic sectors separately. The 
sum of the two is the total ground state energy Eq. Apart from the depen- 
dence on p, Eq is a function of three parameters E = E (A , L, K), where 
we understand our choice of TWF and also optimization of its free param- 
eters. The dependence on K is totally artificial and we are interested in 
E (\ ,L, oo). The dependence on L enters the study of the continuum 
limit of the model. However, here we simply want to test the method and 
work at fixed L = 10. 

For the quadratic case, we show in Fig. (^|) the behavior of E (X , L, K) 
as a function of K at several values of A . For A > —0.5 the K dependence 
is very mild and one can confirm the claim that a positive i?o(^0; L, oo) > 
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Figure 2: Summary plots of the four trial wave function parameters as a function 
of Aq and K. 



is obtained at K — > oo. Instead, for Ao < —0.5, there is a certain 
dependence on K and an extrapolation toward K — > oo must be performed 
to determine the asymptotic limit. Fig. (|j) shows the results of such an 
extrapolation. For Ao > —1.25, we can exclude a zero E (X Q , L, oo); for 
Ao < —1.5 our data are compatible with zero. 

On the other hand, for the cubic potential the dependence of Eq on K 
is relatively mild; in Fig. we show that Eq is compatible with for K > 
500, in full agreement with the expectation of unbroken supersymmetry. 
It should be noticed that bosonic and fermionic contributions to E are of 
the order of 10, and the two are canceling to a precision of 10 -4 . 

The conclusion we can draw from the above numerical data is that 
the presented algorithm behaves in a completely satisfactory way in the 
analysis of the N — 1 Wess-Zumino model, at least for what concerns 
the determination of the ground state energy. In the L — 10 lattice 
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Figure 3: Total energy in log scale as a function of Aq and K. 



quadratic model supersymmetry is broken with no doubts for values of Ao 
down to about A — —1.25. Below this value, E is compatible with zero 
at the statistical level we work. We remark that around A — —0.75 the 
coefficient changes sign, modifying the shape of the TWF and allowing 
for local minima in its (/^-dependent part Sb(^p) at non zero fields; we 
interpret this fact as an interesting signal that some transition is occurring 
and emphasize the important role that trial wave functions play. On the 
other hand, for the cubic model the evidence for unbroken supersymmetry 
is quite convincing. 

Unbroken supersymmetry implies a number of non-trivial Ward identi- 
ties; we monitored several of these, obtaining a pattern of supersymmetry 
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Figure 4: Extrapolation of the total energy in the K — > oo limit. The various 
curves are drawn for several values of Ao- Data points corresponds to K = 100, 
200, 500, 1000, 1700, 3000, 5000, and 10000. 



breaking perfectly consistent with the one obtained from E , although 



with lower numerical precision [11] 



5 Conclusions 

The Hamiltonian approach is a powerful method for Monte Carlo analysis 
of field models. It is well founded and general purpose, but the control 
of systematic errors is fundamental. A good trial wave function can im- 
prove strongly the quality of numerical results and the convergence rate 
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Figure 5: For the cubic potential, extrapolation of the total energy in the 
K -> oo limit. Data points corresponds to K = 100, 200, 500, 1000, and 2000. 



of simulations. The trial wave function is an approximation to the exact 
ground state. As such, it contains important physical informations about 
the model under study. Optimization of (many-parameter) wave functions 
is thus not only motivated by computational needs only. It affords and 
checks analytical insights about the actual ground state. 

Moreover, as we have shown, for certain models in 1 + 1 dimensions, 
the Hamiltonian formalism appears to be the natural framework for lat- 
tice fermions. The treatment of bosons and fermions is nicely symmetric 
and the non-local determinants required in Lagrangian simulations are 
not required at all. The possibility of preserving exactly a 1-dimensional 
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super symmetry algebra is also very appealing. 
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